'''
import sympy as sp
from sympy import numer, denom, Poly
import numpy as np
from numpy import abs, exp, cos, angle
from scipy.signal import residue
from math import factorial
def frac(n,d):
    return sp.Rational(n,d)
S, t = sp.symbols("S t")
'''

'\nimport sympy as sp\nfrom sympy import numer, denom, Poly\nimport numpy as np\nfrom numpy import abs, exp, cos, angle\nfrom scipy.signal import residue\nfrom math import factorial\ndef frac(n,d):\n    return sp.Rational(n,d)\nS, t = sp.symbols("S t")\n'

ECE 2260 Lab 9

Jacob Johnson W01511510

2025 Apr 16

from IPython.display import Math

In this lab, we will use high- and low-pass filters to separate audio frequencies. Circuits like these are commonly used in high-end audio equipment for playing different signals coming from one main signal onto multiple speakers which are tailored to a particular range.

1. Analytic transfer function

We will analytically derive the transfer function $H(S)$ for each circuit. We will be using the sympy python module to handle algebra.

import sympy as sp          # import the module
S, t = sp.symbols("S t")    # create the S and t symbols that we will manipulate

# define helper functions for circuit analysis
def series(*Z):
    return sum(Z)

def parallel(*Z):
    return 1 / sum(1/z for z in Z)

We first implement our components:

L = sp.Rational(1,1_000)            # 1 mH
R = 8                               # 8 Ohm
C_sub = sp.Rational(50,1_000_000)   # 50 uF
C_tweet = sp.Rational(5, 1_000_000) # 5 uF

We can convert these into S-domain:

L = L*S
R = R
C_sub = 1/(C_sub*S)
C_tweet = 1/(C_tweet*S)

Now we can construct our S-domain transfer function:

# subwoofer
H_lp = (parallel(R, C_sub) / series(L, parallel(R, C_sub))).simplify()
display(Math("H_{lp} = " + sp.latex(H_lp)))

# tweeter
H_hp = (parallel(R, L) / series(C_tweet, parallel(R, L))).simplify()
display(Math("H_{hp} = " + sp.latex(H_hp)))

$\displaystyle H_{lp} = \frac{20000000.0}{1.0 S^{2} + 2500.0 S + 20000000.0}$

$\displaystyle H_{hp} = \frac{S^{2}}{S^{2} + 25000.0 S + 200000000}$

Finally, we can take the inverse laplace transform to get the impulse function:

h_lp = sp.inverse_laplace_transform(H_lp, S, t)
display(Math("h_{lp} = " + sp.latex(h_lp)))


h_hp = sp.inverse_laplace_transform(H_hp, S, t)
display(Math("h_{hp} = " + sp.latex(h_hp)))

$\displaystyle h_{lp} = 4657.78007791666 e^{- 1250.0 t} \sin{\left(4293.89100932942 t \right)} \theta\left(t\right)$

$\displaystyle h_{hp} = \left(17008.4012854152 e^{- 12500.0 t} \sin{\left(6614.37827766148 t \right)} - 25000.0 e^{- 12500.0 t} \cos{\left(6614.37827766148 t \right)}\right) \theta\left(t\right) + 1.0 \delta\left(t\right)$

2. Python Simulation

Next, we wil use the scipy.signal python module to generate a Bode plot of the frequency response for each filter.

from matplotlib import pyplot as plt
from scipy import signal
def bode(H, w):
    N_s = np.array(Poly(numer(H), S).all_coeffs()).astype(float)
    D_s = np.array(Poly(denom(H), S).all_coeffs()).astype(float)
    sys = signal.TransferFunction(N_s, D_s)
    w, mag, phase = signal.bode(sys, w)
    return mag, phase

f = np.geomspace(100, 35_000, 100)
w = f*2*np.pi
mag_lp, phase_lp = bode(H_lp, w)
mag_hp, phase_hp = bode(H_hp, w)

fig, axs = plt.subplots(2)
axs[0].plot(f, mag_lp, color="green", label="subwoofer")
axs[0].plot(f, mag_hp, color="blue", label="tweeter")
axs[0].semilogx()
axs[0].set_title("Magnitude of filter transfer function")
axs[0].set_xlabel("Frequency (Hz)")
axs[0].set_ylabel("Gain (dB)")
axs[0].legend()

axs[1].plot(f, phase_lp, color="green", label="subwoofer")
axs[1].plot(f, phase_hp, color="blue", label="tweeter")
axs[1].semilogx()
axs[1].set_title("Phase of filter transfer function")
axs[1].set_xlabel("Frequency (Hz)")
axs[1].set_ylabel("Phase (degrees)")
axs[1].legend()

plt.suptitle("Figure: Bode plot for Subwoofer and Tweeter Networks")
fig.tight_layout()
fig.align_titles()

plt.show()

We can also use these plots to find the cutoff point of our filters:

fc_lp = f[abs(mag_lp - (- 3)).argmin()]
print(f"Cutoff point for subwoofer is {fc_lp} Hz")

fc_hp = f[abs(mag_hp - (- 3)).argmin()]
print(f"Cutoff point for tweeter is {fc_hp} Hz")

Cutoff point for subwoofer is 1066.3659778001622 Hz
Cutoff point for tweeter is 2915.8585332101434 Hz

Now, we reuse code from the previous lab to extract the impulse response of each filter:

# this is my code from last time
def inverse_laplace_lab8( N_s , D_s , t ) :
    pfe = dict()
    K, R, P = residue(N_s, D_s)
    for k, r in zip(K, R):
        if r not in pfe:
            pfe[r] = list()
        pfe[r].append((k.conjugate() if r.imag < 0 else k) / factorial(len(pfe[r])))
    f = np.zeros_like(t)
    for r in pfe:
        for n, k in enumerate(pfe[r]):
            f += (t**n) * abs(k) * exp(r.real*t) * cos(abs(r.imag)*t + angle(k))
    for n, p in enumerate(P[::-1]):
        f += p * t**n
    return f

# this is a wrapper which can accept sympy expressions as the transfer function input
def inverse_laplace(H, t):
    N_s = Poly(numer(H), S).all_coeffs()
    D_s = Poly(denom(H), S).all_coeffs()
    return inverse_laplace_lab8(N_s, D_s, t)

And we plot each response:

duration = .004
sample_rate = 44_100
t = np.linspace(0, duration, int(duration*sample_rate))
h_lp = inverse_laplace(H_lp, t)
h_hp = inverse_laplace(H_hp, t)

plt.plot(t*1000, h_lp, color="green", label="subwoofer")
plt.plot(t*1000, h_hp, color="blue", label="tweeter")
plt.title("Impluse response of filters")
plt.xlabel("t (ms)")
plt.ylabel("Volts")
plt.legend()
plt.show()

In addition to plotting the inpuslse response, we can digitally convolve it with an input signal (Handel's Messiah, for example) to predict what it will sound like when played over the filtered speaker:

import sounddevice as sd
import soundfile as sf
from IPython.display import Audio


filename = "handel.wav"
data, fs = sf.read(filename, dtype = 'float32') # read the file
data = data[:, 0]                               # convert to mono

data_lp = np.convolve(data, h_lp )              # apply low - pass filtering
data_lp /= max(abs(data_lp))                    # normalize, preserve hearing
sf.write("lp.wav", data_lp,  44_100)            # write to .wav file

data_hp = np.convolve(data, h_hp)               # apply high - pass filtering
data_hp /= max ( abs ( data_hp ) )              # normalize, preserve hearing
sf.write("hp.wav", data_hp,  44_100)            # write to .wav file

# play the original
Audio("handel.wav")

# play the low-pass filtered version
Audio("lp.wav")

# play the high-pass filtered version
Audio("hp.wav")

3. LTSpice

We implemented our circuits in LTSpice, then did an AC sweep from $100Hz$ to $35kHz$ to generate a Bode plot for the circuits.

Figure: our filters, modeled in LTSpice.

Figure: the Bode plots for our filters, generated by LTSpice. The dotted lines refer to the phase and use the right y-axis.

import sounddevice as sd
import soundfile as sf

filename = "handel.wav"
data , fs = sf . read ( filename , dtype = 'float32')
data = data [: ,0]
data_lp = np . convolve ( data , h_lp ) # apply low - pass filtering
data_lp /= max ( abs ( data_lp ) ) # normalize , preserve hearing
data_hp = np . convolve ( data , h_hp ) # apply high - pass filtering
data_hp /= max ( abs ( data_hp ) ) # normalize , preserve hearing
# just a break between each of the clips so you can tell
# when there is a transition
sound_break = np . zeros ( fs *2)
# smash everything together into one large array
data_play = np . hstack (( sound_break , data , # original
sound_break , data_lp , # low - pass
sound_break , data , # original
sound_break , data_hp ) ) # high - pass
# play the audio
#sd . play ( data_play , fs )

sd.play(data, fs)

sd.play(data_lp, fs)

sd.play(data_hp, fs)

sd.play(data - data_hp[:len(data)])

length = len(data)//2

data = data[:length]
data_lp = data_lp[:length]
data_hp = data_hp[:length]

sf.write("data.wav", data, 44_100)
sf.write("hp.wav", data_hp,  44_100)
sf.write("lp.wav", data_lp,  44_100)

sys_lp = signal . TransferFunction ([20_000_000],[1, 2_500, 20_000_000])
sys_hp = signal . TransferFunction ([1,0,0],[1, 25_000, 200_000_000])

F = []

for s in [sys_lp, sys_hp]:
    w , mag , phase = signal . bode ( s )
    f0 = (int(w[mag.argmax()] / np.pi / 2))
    fc = (int(w[(mag -- 3).__abs__().argmin()] / np.pi / 2))
    print(fc)
    F.append(f0)
    F.append(fc)
    #print(f0, fc)

print()

F.extend(np.geomspace(100, 35_000, 20).astype(int))

F = sorted(F)

for f in F:
    print(f)

1047
3124

100
136
185
252
343
467
635
642
865
1047
1178
1603
2182
2970
3124
4043
5503
7491
10197
13879
15915
18891
25714
35000

Highpass Hz output phase input 100 .0131 100.5 136 .0165 103 185 .022 111 252 .031 120 343 .0438 123 467 .0643 126 635 .092 124 642 .092 125 865 .122 120 1047 .147 119 1178 .159 116 1603 .207 103 2182 .215 87 2970 .225 72 3124 .227 68 4043 .235 57 5503 .243 44 7491 .249 34 10197 .255 26 13879 .263 20 15915 .269 19 18891 .275 15 25714 .295 12 35000 .324 10

#Highpass
t, mag, phase = np.array([
[100,.0131,100.5],
[136,.0165,103],
[185,.022,111],
[252,.031,120],
[343,.0438,123],
[467,.0643,126],
[635,.092,124],
[642,.092,125],
[865,.122,120],
[1047,.147,119],
[1178,.159,116],
[1603,.207,103],
[2182,.215,87],
[2970,.225,72],
[3124,.227,68],
[4043,.235,57],
[5503,.243,44],
[7491,.249,34],
[10197,.255,26],
[13879,.263,20],
[15915,.269,19],
[18891,.275,15],
[25714,.295,12],
[35000,.324,10]
]).T

db = 20*np.log10(mag)

plt.plot(t, db)
plt.semilogx()
plt.show()

plt.plot(t, phase)
plt.semilogx()
plt.show()

Highpass Hz output phase input 100 .0142 104 2.01 136 .018 106 2.01 185 .024 111 1.97 252 .032 117 1.89 343 .047 123 1.79 467 .067 126 1.63 635 .096 125 1.37 642 .098 127 1.36 865 .132 123 1.11 1047 .154 118 .95 1178 .175 115 .860 1603 .207 101 .660 2182 .235 86 .510 2970 .251 67 .430 3124 .255 64 .410 4043 .263 51 .378 5503 .271 28 .358 7491 .275 29 .350 10197 .275 21 .342 13879 .279 16 .338 15915 .279 13 .334 18891 .279 12 .334 25714 .279 10 .330 35000 .279 6.5 .326

#Highpass],
#[Hz,output,phase,input],

f, output, phase, input = np.array([
[100,.0142,104,2.01],
[136,.018,106,2.01],
[185,.024,111,1.97],
[252,.032,117,1.89],
[343,.047,123,1.79],
[467,.067,126,1.63],
[635,.096,125,1.37],
[642,.098,127,1.36],
[865,.132,123,1.11],
[1047,.154,118,.95],
[1178,.175,115,.860],
[1603,.207,101,.660],
[2182,.235,86,.510],
[2970,.251,67,.430],
[3124,.255,64,.410],
[4043,.263,51,.378],
[5503,.271,28,.358],
[7491,.275,29,.350],
[10197,.275,21,.342],
[13879,.279,16,.338],
[15915,.279,13,.334],
[18891,.279,12,.334],
[25714,.279,10,.330],
[35000,.279,6.5,.326]
]).T

db = 20*np.log10(output/input)

plt.plot(t, db)
plt.semilogx()
plt.show()

plt.plot(t, phase)
plt.semilogx()
plt.show()

Lowpass Hz output phase input 100 .271 -6.7 .330 136 .267 -9 .320 185 .259 -12 .360 252 .243 -17 .285 343 .223 -24 .249 467 .200 -37 .209 635 .167 -59 .171 642 .167 -60 .169 865 .135 -92 .167 1047 .113 -113 .193 1178 .105 -120 .221 1603 .080 -138 .322 2182 .060 -145 .462 2970 .046 -150 .640 3124 .044 -150 .680 4043 .035 -147 .850 5503 .027 -147 1.09 7491 .020 -145 1.33 10197 .016 -132 1.54 13879 .011 -123 1.75 15915 .010 -120 1.81 18891 .010 -120 1.8 25714 .007 1.95 35000 .008 1.98

#Lowpass],
#[Hz,output,phase,input],

f, output, phase, input = np.array([
[100,.271,-6.7,.330],
[136,.267,-9,.320],
[185,.259,-12,.360],
[252,.243,-17,.285],
[343,.223,-24,.249],
[467,.200,-37,.209],
[635,.167,-59,.171],
[642,.167,-60,.169],
[865,.135,-92,.167],
[1047,.113,-113,.193],
[1178,.105,-120,.221],
[1603,.080,-138,.322],
[2182,.060,-145,.462],
[2970,.046,-150,.640],
[3124,.044,-150,.680],
[4043,.035,-147,.850],
[5503,.027,-147,1.09],
[7491,.020,-145,1.33],
[10197,.016,-152,1.54],
[13879,.011,-153,1.75],
[15915,.010,-150,1.81],
[18891,.010,-150,1.8],
[25714,.007,-150,1.95],
[35000,.008,-150, 1.98]]).T

db = 20*np.log10(output/input)

plt.plot(t, db)
plt.semilogx()
plt.show()

plt.plot(t, phase)
plt.semilogx()
plt.show()

460 is lowest I can hear without filter, w LPF, w HPF

13.9k is highest I can hear without filter, 10.2k w LPF, 13.9 w HPF